﻿#pragma once

#include <algorithm>
#include <iostream>
#include <vector>

#include <learn/check.h>
#include <learn/math/math.h>

namespace xi
{

namespace math
{

/**
 * @brief
 *
 * @tparam T 系数类型
 * @tparam Real 参数类型
 */
template <class T, class Real = T> class Polynomial
{
  protected:
    static constexpr auto cmp = Compare{};

    Real m_offset{};
    std::vector<T> m_coefficients;

  public:
    Polynomial(const Polynomial &poly) = default;
    Polynomial(Polynomial &&poly) = default;
    Polynomial &operator=(const Polynomial &poly) = default;
    Polynomial &operator=(Polynomial &&poly) = default;

    /**
     * @brief 限定显式调用，防止 operator* 中的隐式转换
     *
     * @param[in] deg
     */
    explicit Polynomial(std::size_t deg) : m_coefficients(deg + 1, T{})
    {
    }

    /**
     * @brief 系数构造函数
     *
     * @param[in] coefficients
     */
    explicit Polynomial(std::vector<T> &&coefficients) : m_coefficients(std::move(coefficients))
    {
    }

    /**
     * @brief 比较
     *
     * @param[in] poly
     * @return true
     * @return false
     */
    bool operator==(const Polynomial &poly) const
    {
        if (degree() != poly.degree())
            return false;

        for (std::size_t i = 0; i < degree() + 1; i++)
            if (!cmp.test<Real>(m_coefficients[i], poly.m_coefficients[i]))
                return false;

        return true;
    }

    /**
     * @brief 参数偏移
     *
     * @param[in] offset
     * @return Polynomial
     */
    Polynomial &operator>>(Real offset)
    {
        m_offset += offset;
        return *this;
    }

    /**
     * @brief 获得偏移值
     *
     * @return Real
     */
    Real offset() const
    {
        return m_offset;
    }

    /**
     * @brief 获得次数
     *
     * @return std::size_t
     */
    constexpr std::size_t degree() const noexcept
    {
        return m_coefficients.size() - 1;
    }

    /**
     * @brief 升次
     *
     * @param[in] d
     * @return true
     * @return false
     */
    bool elevate(std::size_t d)
    {
        if (d <= degree())
            return false;

        // 高位填 0
        m_coefficients.resize(d + 1, T{});
        return true;
    }

    /**
     * @brief 多项式降阶
     *
     * @param[in] d
     * @return true
     * @return false
     */
    bool reduce(std::size_t d)
    {
        if (d >= degree())
            return false;

        // 低位截断
        m_coefficients.resize(d + 1, T{});
        return true;
    }

    /**
     * @brief 获得系数
     *
     * @return const std::vector<T>&
     */
    const std::vector<T> &coefficients() const noexcept
    {
        return m_coefficients;
    }

    /**
     * @brief 引用系数
     *
     * @return std::vector<T>&
     */
    std::vector<T> &coefficients() noexcept
    {
        return m_coefficients;
    }

    /**
     * @brief 获得系数
     *
     * @param[in] i
     * @return T
     */
    T coefficient(std::size_t i) const
    {
        return m_coefficients[i];
    }

    /**
     * @brief 引用系数
     *
     * @param[in] i
     * @return T&
     */
    T &coefficient(std::size_t i)
    {
        return m_coefficients[i];
    }

#define POLY_OP(func, op)                                                                                              \
     Polynomial func(const Polynomial &poly) const                                                                \
    {                                                                                                                  \
        Polynomial res(std::max(degree(), poly.degree()));                                                             \
        auto &coeff = res.coefficients();                                                                              \
        auto &coeff1 = m_coefficients;                                                                                 \
        auto &coeff2 = poly.coefficients();                                                                            \
        for (std::size_t i = 0; i < coeff1.size(); i++)                                                                \
            coeff[i] = coeff1[i];                                                                                      \
        for (std::size_t i = 0; i < coeff2.size(); i++)                                                                \
            coeff[i] = coeff[i] op coeff2[i];                                                                          \
        return res;                                                                                                    \
    }

    POLY_OP(operator+, +)
    POLY_OP(operator-, -)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
    const Polynomial &func(const Polynomial &poly)                                                                     \
    {                                                                                                                  \
        elevate(poly.degree());                                                                                        \
        auto &coeff = m_coefficients;                                                                                  \
        auto &coeff2 = poly.coefficients();                                                                            \
        for (std::size_t i = 0; i < coeff2.size(); i++)                                                                \
            coeff[i] op coeff2[i];                                                                                     \
        return *this;                                                                                                  \
    }

    POLY_OP(operator+=, +=)
    POLY_OP(operator-=, -=)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
     Polynomial func(T k) const                                                                                   \
    {                                                                                                                  \
        Polynomial res(degree());                                                                                      \
        for (std::size_t i = 0; i < m_coefficients.size(); i++)                                                        \
            res.m_coefficients[i] = m_coefficients[i] op k;                                                            \
        return res;                                                                                                    \
    }

    POLY_OP(operator*, *)
    POLY_OP(operator/, /)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
    const Polynomial &func(T k)                                                                                        \
    {                                                                                                                  \
        for (std::size_t i = 0; i < m_coefficients.size(); i++)                                                        \
            m_coefficients[i] op k;                                                                                    \
        return *this;                                                                                                  \
    }

    POLY_OP(operator*=, *=)
    POLY_OP(operator/=, /=)
#undef POLY_OP

    /**
     * @brief 多项式乘法
     *
     * @param[in] poly
     * @return Polynomial
     */
     Polynomial operator*(const Polynomial &poly) const
    {
        Polynomial res(degree() + poly.degree());
        for (std::size_t i = 0; i < degree() + 1; i++)
            for (std::size_t j = 0; j < poly.degree() + 1; j++)
                res.m_coefficients[i + j] += m_coefficients[i] * poly.m_coefficients[j];
        return res;
    }

    /**
     * @brief 多项式乘法
     *
     * @param[in] poly
     * @return Polynomial&
     */
    const Polynomial &operator*=(const Polynomial &poly)
    {
        *this = *this * poly;
        return *this;
    }

    /**
     * @brief 多项式除法
     *
     * @param[in] poly
     * @return Polynomial
     */
     Polynomial operator/(const Polynomial &poly) const
    {
        // 次数不够，返回 0 多项式（这里是唯一允许次数为负的情况）
        if (degree() < poly.degree())
            return Polynomial(-1);

        Polynomial res(degree() - poly.degree());
        res.m_coefficients[degree() - poly.degree()] = m_coefficients.back() / poly.m_coefficients.back();

        // 获得余式
        Polynomial tmp = *this - res * poly;

        // 减去最高项后需要降阶，然后递归计算商
        tmp.reduce(degree() - 1);
        tmp = tmp / poly;

        // 商式赋值
        for (std::size_t i = 0; i < tmp.degree() + 1; ++i)
            res.m_coefficients[i] = tmp.m_coefficients[i];

        return res;
    }

    /**
     * @brief 多项式除法
     *
     * @param[in] poly
     * @return Polynomial&
     */
    const Polynomial &operator/=(const Polynomial &poly)
    {
        *this = *this / poly;
        return *this;
    }

    /**
     * @brief 求值
     *
     * @param[in] x
     * @return T
     */
    T value(Real x) const
    {
        std::size_t O = degree();
        if (O < 0)
            return T{};

        // 获得最高次系数，使用累乘方法提高效率
        T y = m_coefficients[O];
        x -= m_offset;
        for (int i = O - 1; i >= 0; i--)
            y = y * x + m_coefficients[i];

        return y;
    }

    /**
     * @brief 多项式定积分
     *
     * @param[in] a
     * @param[in] b
     * @return T
     */
    T integral(Real a, Real b) const
    {
        T I = T{};
        a -= m_offset;
        b -= m_offset;
        for (std::size_t i = 0; i < degree() + 1; i++)
            I += m_coefficients[i] * (std::pow(b, i + 1) - std::pow(a, i + 1)) / (i + 1);
        return I;
    }

    /**
     * @brief 多项式不定积分
     *
     * @return Polynomial
     */
    Polynomial integral() const
    {
        Polynomial poly(degree() + 1);
        for (std::size_t i = 0; i < degree() + 1; i++)
            poly.coefficients()[i + 1] = m_coefficients[i] / (i + 1);
        return poly >> m_offset;
    }

    /**
     * @brief 多项式求导
     *
     * @return Polynomial
     */
    Polynomial derivative() const
    {
        if (degree() == 0)
            return Polynomial(0);

        Polynomial poly(degree() - 1);
        for (std::size_t i = 1; i < degree() + 1; i++)
            poly.coefficients()[i - 1] = m_coefficients[i] * i;
        return poly >> m_offset;
    }

    /**
     * @brief 产生节点多项式
     *
     * @param[in] c
     * @return Polynomial
     */
    static Polynomial nodes(const std::vector<Real> &c)
    {
        // 直接写 {1} 似乎会被错误识别
        Polynomial poly(0);
        poly.coefficient(0) = T{1};
        for (std::size_t i = 0; i < c.size(); i++)
            poly = poly * Polynomial({T{-c[i]}, T{1}});

        return poly;
    }

    /**
     * @brief 已知节点 c0 c1 ... cn ，获得 cj 对应的初等多项式
     *
     * @param[in] c
     * @param[in] j
     * @return Polynomial
     */
    static Polynomial lagrange(const std::vector<Real> &c, std::size_t j)
    {
        // 直接写 {1} 似乎会被错误识别
        Polynomial poly(0);
        poly.coefficient(0) = T{1};
        for (std::size_t k = 0; k < c.size(); k++)
            if (k != j)
                poly = poly * Polynomial({T{-c[k]}, T{1}}) / (c[j] - c[k]);

        return poly;
    }

    /**
     * @brief (-1,1) 上的 Gauss-Legendre 多项式，用于数值积分
     * 三项递推式 i * p_{i-1} - (2 * i + 1) * x * p_i + (i + 1) * p_{i+1} = 0
     *
     * @param[in] s
     * @return Polynomial
     */
    static Polynomial legendre(std::size_t s)
    {
        // 直接写 {1} 似乎会被错误识别
        Polynomial f0(0);
        f0.coefficient(0) = T{1};
        if (s > 0)
        {
            Polynomial x({T{0}, T{1}});
            Polynomial f1 = x, f2 = x;
            for (std::size_t i = 1; i < s; i++)
            {
                f2 = x * f1 * (2 * i + 1) / (i + 1) - f0 * i / (i + 1);
                f0 = f1;
                f1 = f2;
            }
            f0 = f2;
        }
        return f0;
    }

    /**
     * @brief (0,1) 上的首一平移 Gauss-Legendre 多项式，用于 Runge-Kutta 积分
     *
     * @param[in] s
     * @return Polynomial
     */
    static Polynomial monic_legendre(std::size_t s)
    {
        // 注意不要使用 0 多项式
        Polynomial poly(0);
        for (std::size_t j = 0; j <= s; j++)
        {
            Polynomial v(j);
            v.coefficient(j) = T{Real(((s - j) % 2 == 0 ? 1 : -1))} * binomial(s, j) * binomial(s + j, j);
            poly = poly + v;
        }
        poly *= std::pow(factorial(s), 2) / factorial(2 * s);
        return poly;
    }

    /**
     * @brief 输出函数
     *
     * @param[in] os
     * @param[in] poly
     * @return std::ostream&
     */
    friend std::ostream &operator<<(std::ostream &os, const Polynomial &poly)
    {
        auto &coeff = poly.coefficients();

        os << "f(x) = ";
        for (std::size_t i = 0; i < coeff.size(); i++)
        {
            if (i == 0)
                os << coeff[i];
            else
                os << " + " << coeff[i] << "x^" << i;
        }
        return os;
    }
};

template <class T, class Real = T> class PiecewisePolynomial
{
  protected:
    static constexpr auto cmp = Compare{};

    std::vector<Real> m_interval;            // 分段区间端点，区间左闭右开
    std::vector<Polynomial<T, Real>> m_poly; // 多项式列表，数量是端点数 +1，包括左右半开区间上的多项式

  public:
    PiecewisePolynomial(const PiecewisePolynomial &) = default;
    PiecewisePolynomial(PiecewisePolynomial &&) = default;
    PiecewisePolynomial &operator=(const PiecewisePolynomial &) = default;
    PiecewisePolynomial &operator=(PiecewisePolynomial &&) = default;

    /**
     * @brief 空构造
     *
     */
    explicit PiecewisePolynomial() : m_poly({Polynomial<T, Real>(0)})
    {
    }

    /**
     * @brief 单多项式
     *
     * @param[in] poly
     */
    explicit PiecewisePolynomial(Polynomial<T, Real> &&poly) : m_poly({std::move(poly)})
    {
    }

    /**
     * @brief 只设置区间，每段区间上是零
     *
     * @param[in] interval
     */
    explicit PiecewisePolynomial(std::vector<Real> &&interval)
        : m_interval(std::move(interval)), m_poly(m_interval.size() + 1, Polynomial<T, Real>(0))
    {
    }

    /**
     * @brief 设置区间和多项式
     *
     * @param[in] interval
     * @param[in] poly
     */
    PiecewisePolynomial(std::vector<Real> &&interval, std::vector<Polynomial<T, Real>> &&poly)
        : m_interval(std::move(interval)), m_poly(std::move(poly))
    {
        // 检查赋值是否有效
        ASSERT_INFO(m_interval.size() + 1 == m_poly.size(), "PiecewisePolynomial: invalid input size");
    }

    /**
     * @brief 比较
     *
     * @param[in] poly
     * @return true
     * @return false
     */
    bool operator==(const PiecewisePolynomial &poly) const
    {
        return m_interval == poly.m_interval && m_poly == poly.m_poly;
    }

#define POLY_OP(func, op)                                                                                              \
     PiecewisePolynomial func(const PiecewisePolynomial &poly) const                                              \
    {                                                                                                                  \
        ASSERT_INFO(m_interval.size() == poly.m_interval.size(), "PiecewisePolynomial: invalid input size");           \
        for (std::size_t i = 0; i < m_interval.size(); ++i)                                                            \
            ASSERT_INFO(std::abs(m_interval[i] - poly.m_interval[i]) < std::numeric_limits<Real>::epsilon(),           \
                        "PiecewisePolynomial: invalid input interval");                                                \
        auto interval = m_interval;                                                                                    \
        PiecewisePolynomial res(std::move(interval));                                                                  \
        for (std::size_t i = 0; i < m_poly.size(); ++i)                                                                \
            res.m_poly[i] = m_poly[i] op poly.m_poly[i];                                                               \
        return res;                                                                                                    \
    }

    POLY_OP(operator+, +)
    POLY_OP(operator-, -)
    POLY_OP(operator*, *)
    POLY_OP(operator/, /)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
    const PiecewisePolynomial &func(const PiecewisePolynomial &poly)                                                   \
    {                                                                                                                  \
        ASSERT_INFO(m_interval.size() == poly.m_interval.size(), "PiecewisePolynomial: invalid input size");           \
        for (std::size_t i = 0; i < m_interval.size(); ++i)                                                            \
            ASSERT_INFO(std::abs(m_interval[i] - poly.m_interval[i]) < std::numeric_limits<Real>::epsilon(),           \
                        "PiecewisePolynomial: invalid input interval");                                                \
        for (std::size_t i = 0; i < m_poly.size(); ++i)                                                                \
            m_poly[i] op poly.m_poly[i];                                                                               \
        return *this;                                                                                                  \
    }

    POLY_OP(operator+=, +=)
    POLY_OP(operator-=, -=)
    POLY_OP(operator*=, *=)
    POLY_OP(operator/=, /=)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
     PiecewisePolynomial func(T k) const                                                                          \
    {                                                                                                                  \
        auto interval = m_interval;                                                                                    \
        PiecewisePolynomial res(std::move(interval));                                                                  \
        for (std::size_t i = 0; i < m_poly.size(); ++i)                                                                \
            res.m_poly[i] = m_poly[i] op k;                                                                            \
        return res;                                                                                                    \
    }

    POLY_OP(operator*, *)
    POLY_OP(operator/, /)
#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
    const PiecewisePolynomial &func(T k)                                                                               \
    {                                                                                                                  \
        for (auto &poly : m_poly)                                                                                      \
            poly op k;                                                                                                 \
        return *this;                                                                                                  \
    }

    POLY_OP(operator*=, *=)
    POLY_OP(operator/=, /=)
#undef POLY_OP

    /**
     * @brief 获得多项式数量
     *
     * @return std::size_t
     */
    constexpr std::size_t poly_num() const noexcept
    {
        return m_poly.size();
    }

    /**
     * @brief 获得端点
     *
     * @return const std::vector<Real>&
     */
    const std::vector<Real> &intervals() const noexcept
    {
        return m_interval;
    }

    /**
     * @brief 获得多项式
     *
     * @return const std::vector<Polynomial<T, Real>>&
     */
    const std::vector<Polynomial<T, Real>> &polynomials() const noexcept
    {
        return m_poly;
    }

    /**
     * @brief 获得多项式
     *
     * @return std::vector<Polynomial<T, Real>>&
     */
    std::vector<Polynomial<T, Real>> &polynomials() noexcept
    {
        return m_poly;
    }

    /**
     * @brief 获得指定索引的次数
     *
     * @param[in] i
     * @return std::size_t
     */
    constexpr std::size_t degree(std::size_t i) const noexcept
    {
        ASSERT_INFO(i < m_poly.size(), "PiecewisePolynomial: invalid index");
        return m_poly[i].degree();
    }

    /**
     * @brief 获得指定索引的多项式
     *
     * @param[in] i
     * @return const Polynomial&
     */
    const Polynomial<T, Real> &operator[](std::size_t i) const
    {
        ASSERT_INFO(i < m_poly.size(), "PiecewisePolynomial: invalid index");
        return m_poly[i];
    }

    /**
     * @brief 获得指定索引的多项式引用
     *
     * @param[in] i
     * @return Polynomial&
     */
    Polynomial<T, Real> &operator[](std::size_t i)
    {
        ASSERT_INFO(i < m_poly.size(), "PiecewisePolynomial: invalid index");
        return m_poly[i];
    }

    /**
     * @brief 获得传入点对应的多项式（左闭右开区间）
     *
     * @param[in] t
     * @return const Polynomial&
     */
    const Polynomial<T, Real> &polynomial(Real t) const
    {
        if (m_interval.empty())
            return m_poly.front();

        // 先排除超出区间范围的部分
        if (t < m_interval.front())
            return m_poly.front();
        if (t >= m_interval.back())
            return m_poly.back();
        if (t == m_interval.front())
            return m_poly[1];

        auto it = std::lower_bound(m_interval.begin(), m_interval.end(), t);
        return m_poly[it - m_interval.begin()];
    }

    /**
     * @brief 多项式求值
     *
     * @param[in] x
     * @return T
     */
    T value(Real x) const
    {
        return polynomial(x).value(x);
    }

    /**
     * @brief 在右端添加端点，增加一个零多项式
     *
     * @param[in] t
     * @return Polynomial<T, Real>&
     */
    Polynomial<T, Real> &push_back(Real t)
    {
        // 如果区间端点不是空，就检查是否大于最大的端点，要确保不能重合
        ASSERT_INFO(m_interval.empty() || t > m_interval.back() + std::numeric_limits<Real>::epsilon(),
                    "PiecewisePolynomial: invalid input node");

        m_interval.emplace_back(t);
        m_poly.emplace_back(0);
        return m_poly.back();
    }

    /**
     * @brief 在左端添加端点，增加一个零多项式
     *
     * @param[in] t
     * @return Polynomial<T, Real>&
     */
    Polynomial<T, Real> &push_front(Real t)
    {
        ASSERT_INFO(m_interval.empty() || t < m_interval.front() - std::numeric_limits<Real>::epsilon(),
                    "PiecewisePolynomial: invalid input node");

        m_interval.emplace(m_interval.begin(), t);
        m_poly.emplace(m_poly.begin(), 0);
        return m_poly.front();
    }

    /**
     * @brief 求导
     *
     * @return PiecewisePolynomial
     */
    PiecewisePolynomial derivative() const
    {
        auto interval = m_interval;
        PiecewisePolynomial res(std::move(interval));
        for (std::size_t i = 0; i < m_poly.size(); ++i)
            res.m_poly[i] = m_poly[i].derivative();
        return res;
    }

    /**
     * @brief 不定积分
     *
     * @return PiecewisePolynomial
     */
    PiecewisePolynomial integral() const
    {
        auto interval = m_interval;
        PiecewisePolynomial res(std::move(interval));
        for (std::size_t i = 0; i < m_poly.size(); ++i)
            res.m_poly[i] = m_poly[i].integral();
        return res;
    }

    /**
     * @brief 定积分
     *
     * @param[in] a
     * @param[in] b
     * @return T
     */
    T integral(Real a, Real b) const
    {
        ASSERT_INFO(a < b, "PiecewisePolynomial: invalid input interval");

        int begin = std::lower_bound(m_interval.begin(), m_interval.end(), a) - m_interval.begin();
        int end = std::lower_bound(m_interval.begin(), m_interval.end(), b) - m_interval.begin();

        if (begin == end)
            return m_poly[begin].integral(a, b);

        T res = m_poly[begin].integral(a, m_interval[begin]);
        for (int i = begin + 1; i < end; i++)
            res += m_poly[i].integral(m_interval[i - 1], m_interval[i]);
        res += m_poly[end].integral(m_interval[end - 1], b);
        return res;
    }

    /**
     * @brief 检查端点处的连续性
     *
     * @param[in] i
     * @return std::size_t
     */
    std::size_t continuity(std::size_t i) const
    {
        ASSERT_INFO(i < m_interval.size(), "PiecewisePolynomial: invalid index");
        if (poly_num() == 1)
            return std::numeric_limits<std::size_t>::infinity();

        // 左右两个多项式
        Polynomial left = m_poly[i];
        Polynomial right = m_poly[i + 1];

        std::size_t min = std::min(left.degree(), right.degree());
        Real t = m_interval[i];

        // 标记连续性
        std::size_t c = -1;
        for (std::size_t j = 0; j < min; j++)
        {
            if (!cmp.test<Real>(left.value(t), right.value(t)))
                break;
            left = left.derivative();
            right = right.derivative();
            c++;
        }
        return c;
    }

    /**
     * @brief 整体连续性
     *
     * @return std::size_t
     */
    std::size_t continuity() const
    {
        if (poly_num() == 1)
            return std::numeric_limits<std::size_t>::infinity();

        std::size_t c = continuity(0);
        for (std::size_t i = 1; i < m_interval.size(); i++)
            c = std::min(c, continuity(i));
        return c;
    }

    /**
     * @brief 流式输出
     *
     * @param[in] os
     * @param[in] ppoly
     * @return std::ostream&
     */
    friend std::ostream &operator<<(std::ostream &os, const PiecewisePolynomial &ppoly)
    {
        auto interval = ppoly.m_interval;
        interval.emplace(interval.begin(), std::numeric_limits<Real>::lowest());
        interval.push_back(std::numeric_limits<Real>::infinity());

        os << "f(x) = {" << std::endl;
        for (std::size_t i = 0; i < ppoly.m_poly.size(); ++i)
        {
            os << "f_" << i << "(x) = " << ppoly.m_poly[i] << " in [" << interval[i] << ", " << interval[i + 1] << "]"
               << std::endl;
        }
        os << "}";
        return os;
    }
};

/**
 * @brief 合并两个 PiecewisePolynomial
 *
 * @tparam T
 * @tparam Real
 * @param[in] p1
 * @param[in] p2
 * @return PiecewisePolynomial<T, Real>
 */
template <class T, class Real = T>
inline PiecewisePolynomial<T, Real> merge(const PiecewisePolynomial<T, Real> &p1,
                                          const PiecewisePolynomial<T, Real> &p2)
{
    if (p1.poly_num() == 1 && p2.poly_num() == 1)
        return PiecewisePolynomial<T, Real>(p1[0] + p2[0]);

    // 空多项式
    PiecewisePolynomial<T, Real> res;

    // 储存所有端点，然后排序
    std::vector<Real> interval;
    interval.insert(interval.end(), p1.intervals().begin(), p1.intervals().end());
    interval.insert(interval.end(), p2.intervals().begin(), p2.intervals().end());
    std::sort(interval.begin(), interval.end());

    // 推入这些端点
    res.push_back(interval.front());
    for (std::size_t i = 1; i < interval.size(); i++)
    {
        if (std::abs(interval[i] - interval[i - 1]) < std::numeric_limits<Real>::epsilon())
            continue;
        res.push_back(interval[i]);
    }

    // 获取对应的多项式相加
    interval.push_back(interval.back() + 1);
    interval.emplace(interval.begin(), interval.front() - 1);
    for (std::size_t i = 0; i < res.poly_num(); i++)
    {
        Real mid = (interval[i] + interval[i + 1]) / 2;
        res[i] = p1.polynomial(mid) + p2.polynomial(mid);
    }
    return res;
}

} // namespace math

} // namespace xi
